Phonon mediated tunneling into graphene 
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Recent scanning tunneling spectroscopy experiments[l, 2] on graphene reported an unexpected 
gap of about ±60 meV around the Fermi level. Here, we give a theoretical investigation explaining 
the experimentally observed spectra and confirming the phonon mediated tunneling as the reason for 
the gap: We study the real space properties of the wave functions involved in the tunneling process 
by means of ab-initio theory and present a model for the electron-phonon interaction, which couples 
the graphene's Dirac electrons with quasi free electron states at the Brillouin zone center. The 
self-energy associated with this electron-phonon interaction is calculated and its effects on tunneling 
into graphene are discussed. In particular, good agreement of the tunneling density of states within 
our model and the experimental[l, 2] AI/AU spectra is found. 



Graphene, the two dimensional allotrope of carbon [3], 
is famous for its electrons having vanishing effective mass 
[4, 5]. At two non-equivalent points at the corner of the 
Brillouin zone, K and K', the linearly dispersing valence 
and conduction band touch. So, electrons in graphene be- 
have like massless Dirac fermions with the speed of light 
being replaced by the Fermi velocity x>f ~ c/300. This 
yields an interesting analogy between condensed matter 
and high energy physics [6], which has been extensively 
studied during the last years [7]. We argue, however, that 
there is a decisive difference between the Dirac fermions 
in graphene and Dirac particles studied in the context 
of high energy physics: The electrons in graphene are 
embedded to a real material which interacts with these 
electrons. This "real material background" has proven to 
introduce spatial charge inhomogeneities due to the im- 
perfections present in any solid [8] : The imperfections in- 
clude charged impurities or structural corrugations. Nei- 
ther free standing graphene nor graphene on a substrate 
is perfectly flat [9, 10], which leads to effective gauge 
fields acting on the Dirac electrons [11]. 

Recently reported scanning tunneling spectroscopy 
(STS) experiments on graphene stress this real material 
backround even more [1, 2], It is not the Dirac point 
with its linearly vanishing density of density of states 
which causes the most prominent feature in the mea- 
sured dl/dll spectra, but a gap of ±60meV pinned to 
the Fermi level [2]. After initial speculations on the na- 
ture of this gap like feature including substrate and elec- 
tric field effects [1], the most recent experiments indicate 
that this gap is caused by the opening of an inelastic tun- 
neling channel due to graphene's out-of-plane at K and 
K' [2]. 

Here, we give a detailled theoretical study of the latter 
scenario. Firstly, we investigate the effect of the K/K' 
out-of-plane phonons on the electronic wave functions 
and their decay in the vacuum by means of first principles 
theory. We show why the electron-phonon coupling has a 
large impact on STS as seen in Refs. [1, 2]. Motivated by 



this insight, we present a simple model of graphene's elec- 
trons being coupled to the out-of-plane phonons. Within 
this model the electron self-energy, the total (DOS) and 
the tunneling density of states (TDOS) are calculated. 
It turns out that the total DOS and the TDOS, which 
is significant for STS, differ strongly: At low energies, 
the total DOS is dominated by the T^-shape from the 
Dirac electron's DOS. However, the TDOS recovers the 
gapped experimental d//d£7 spectra. Here, the inelastic 
channel is strongly enhanced in the TDOS due to a very 
general mechanism of band mixing. This mechanism is 
not limited to dynamic processes but is also expected to 
occur, e.g., near short range corrugations of the graphene 
lattice. 

To address the effects of the K/K' phonons on the elec- 
tronic wave functions of graphene, we performed density 
functional theory calculations within the framework of 
the local density approximation using the Vienna Ab Ini- 
tio Simulation Package (VASP) [12] with the projector 
augmented wave (PAW) [13, 14] basis sets. The corre- 
sponding plane wave expansions were cut-off at 928 eV 
and the Brillouin zone integrations were carried out with 
the tetrahedron-method on k-meshes denser than 20 x 20 
when folded back to the simple graphene unit cell. The 
vertical extension of this cell was chosen to be 24A. 

The band structure of flat single layer graphene ob- 
tained in this way is shown in the upper part of Fig. 1 . 
The 7r-bands (fat red dots) exhibit the well known shape 
intersecting the Fermi level at the K point with linear 
dispersion in the vicinity. At the zone center T, a quasi 
continuum of free electron like bands starts 3.3 eV above 
the Fermi level. These bands correspond to the interlayer 
band of graphite, which plays a crucial role in graphite in- 
tercalation compounds [15]. In contrast to graphite, the 
quasi free electron bands, here, are not limited by other 
graphene layers on top and extend far into the vacuum 
above the sheet. 

To quantify the decay of the wave functions of dif- 
ferent bands in the vacuum, we calculated the laterally 
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Figure 1: (Color online) Upper panel: Band structure of 
graphene. The 7r-bands inhibiting the Dirac electrons are 
marked with fat red dots. Instead of the graphite interlayer 
band, in graphene a quasi-continuum of nearly free electron 
states begins 3.3 eV above the Fermi level at the F-point. 
Lower panel: Decay of the electronic wave functions ^(z) as 
a function of hight z above the graphene sheet. The laterally 
averaged density |^(z)| 2 is shown for wave functions near the 
Dirac point of flat graphene (solid line) and in presence of the 
frozen phonon (dashed line). The density belonging to the 
lowest quasi free electron band at F printed as dotted line. 



averaged charge density |^(z)| 2 of these wave functions 
as a function of distance z from the sheet. As shown 
in Fig. 1, lower panel, the wave functions at the Dirac 
point of flat graphene have their maximum at z = 0.4 A 
above the sheet. For z > 1.5A, their probability den- 
sity starts decaying into the vacuum dominated by an 
exponential with decay constant A^- 1 = 3.9A -1 , which 
has been obtained from a linear fit of log ^^(z)! 2 in 
this region. This decay constant follows from the sin- 
gle particle Schrodinger equation in presence of a con- 
stant potential given by the graphene work function [16]: 

A~ 4 = fc|| + 2rrif(j)/h 2 , where the in- plane crystal mo- 
mentum is |fcn| = 1.7A -1 for the electrons at K/K' and 
we obtained <fi = 4.5 eV for the work function from LDA. 

States at the T point of the Brillouin zone, fc| = 0, 
decay much slower into the vacuum: The lowest nearly 
free electron state is extremely long ranged (see Fig 1., 
lower panel) and spreads out into the vacuum with decay 
constant A -1 = 1.0A _1 . Near the Fermi level Ef, slowly 
decaying electronic states are generated by out-of-plane 
K- and K'-phonons, which mix the nearly free electron 
bands at T with the Dirac-like w bands at K. To illus- 
trate this point, consider a frozen phonon with amplitude 



Az = O.lA and a wave vector connecting K/K'. In pres- 
ence of this phonon the Dirac bands formerly at K and 
K' are folded back to the T point of the x % /3) R30° 
super cell inhibiting the frozen phonon. One out of the 
Dirac bands touching the Fermi level at T hybridizes with 
the long range quasi free electron states upon formation 
of the frozen phonon. This results in the long range decay 
of the corresponding charge density, |\E'r(-2 ; )| 2 ~ e~ A r z 
with Ap 1 = 2.2A -1 , being approximately twice slower 
than for the unperturbed Dirac bands. 

Experimentally, the decay of the wave functions in- 
volved in the tunneling processes has been obtained from 
the dependence of the tunneling current on the tip- 
sample distance [2]. Different decay constants Ajl ~ 
4A" 1 and A^t ~ 2A- 1 were measured inside and out- 
side the gap [2], respectively, which nicely coincide with 
the theoretical values A^- 1 and Ap . 

According to Tersoff-Hamann theory, tunneling cur- 
rents are determined by the tunneling density of states, 
which is the LDOS at the location of the tip - usu- 
ally on the order 5A above the sample. Thus, for each 
band, the DOS has to be weighted with squared am- 
plitude of the corresponding wave function about 5A 
above the sample - which can be approximately writ- 
ten as dl/dU ~ \^ r \ 2 N r {E) + \^ K \ 2 N K {E). With the 
decay constants from above for the Dirac bands of flat 
graphene and in presence of the phonon, the ratio of the 
exponential prefactors entering the tunneling matrix el- 
ements is l^r/^ifl 2 oc e 17A z . At z = 5A this ratio 
is on the order of 10 4 . Thus, the effect of mixing of the 
free electron bands and the 7r-bands by electron-phonon 
coupling on STS spectra will be greatly enhanced by this 
factor of 10 4 . 

In the following, we present a simple model to under- 
stand how electron-phonon coupling due to the K/K' out 
of plane phonons affects STS: As the previous discussion 
showed, most important for these experiment will be the 
mixing of the quasi free electron bands at T with the 
Dirac fermion bands at K/K'. The former bands have 
their minimum E a = 3.3 eV above the Fermi level and 
may be approximated by a flat band H a = Y] q S q E a d q 
in the vicinity of T. Here, d q is the annihilation op- 
erator of an electron with crystal momentum q in the 
nearly free electron band. The two tt bands give rise 
to graphene's linear density of states close to the Dirac 
point and their Hamiltonian may be written as = 
Si/=± q vt { ( i)Aj,q c 'j,q-i where the index v = +(— ) denotes 
the conduction (valence) band and c u . q is the annihilator 
of an electron in this 7r-band with momentum q. Close 
to the Dirac points K and K', i.e. q = ±K + k with 
|fc| <C \K\, the dispersion is linear e(k) « ?iuf|fc|. 

Out-of-plane phonons scatter electrons between the it 
and quasi free electron bands with the electron-phonon- 
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interaction reading as 

V = X J2 (4+,^,9 + 4,fc+,^)K + a -fe)' (!) 

where annihilates an out-of-plane phonon carrying 
crystal momentum k. As stated above, the phonon 
modes at K and K' will be the most important con- 
tributors to inelastic tunneling signals. Around these 
points, their Hamiltonian can be approximated by -ff p h = 
hu! a J2k a k a k WR *h = 67meV being the energy of the 
out-of-plane phonons at K and K'. [17] 

In this model dealing with three electronic bands, the 
non-interacting electron Green function is a diagonal 3x3 
matrix and reads as 



G°(q, iuj n ) 



( 



V 



iu> n —E a 










\ 



(2) 



where uj n are fermionic Matsubara frequencies. With 



the non-interacting phonon Green function D°(k, iSl m ) = 
D°(iQ m ) = — q:?""^ an d the electron-phonon interac- 
tion from Eqn. (1) transformed to the matrix form of 

/ A \ /OOA 
Eqn. (2), M_ = A and M+ = 

\0 00/ \ A 

the electronic self-energy is 



E(iw„) = ~ D(iSl m )M v G°(r = 0,iu n -iQ m )M v , 

(3) 

where G°(r — 0, iu n ) is the non-interacting local Green's 
function and the independence of D° and M„ of k has 
been exploited. This self-energy is diagonal and, for 
energies u <C W small as compared to the Dirac elec- 
tron bandwidth W « 6 eV, its components read as 
Ei,i(w + iS) = Si ^w + iS) + iE'/^w + iS), where 
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and 



2ttA 



Tt'l^u + iS) = - LUa)\u + M~ sign(w)w a | 



The 7T block is 



E^w + iJ) = A 2 



11- E a 



(5) 



x 2 



with i = 1,2 and causes a rigid shift of the band energies 
in regime important for STS. 

With the interacting Green function being G _1 (p, w) = 

G^ -* (p, w) — E(o>) momentum space integration yields 
the total density of states N(u) = — ^ Tr ImG(r = 
0,w + For fixed chemical potential /i = —OAeV, 
the upper panel in Fig. 2 shows how the total DOS is 
modified to the electron-phonon interaction as a func- 
tion of the coupling strength A. The total DOS and the 
DOS of the 7r-bands (not shown here) are virtually indis- 
tinguishable and exhibit the " V" -shape characteristic for 
Dirac fermions. The main modification of total DOS is 
a shift to slightly lower energies with increasing coupling 
strength due to the 7r-block of T,(uj + i5) — a consequence 
of level repulsion of the n and the nearly free electron 
bands, well known to occur in second order perturbation 
theory. 



However, the spectral properties of the nearly free 
electron channel are strongly altered at low energies 
\u\ <C W. The density of states in this channel reads 
as 



N r (oj) = - 



it\uj — E a — Ei,i(cj + id)\ 2 
~ 6(M - uj a ) \u + fi - sign(w)w a | . (7) 

Without electron-phonon interaction, the DOS of the 
nearly free electron states vanishes for ui < E a . But, as 
soon as this interaction becomes effective and the elec- 
tron's energy u> > u> a exceeds the phonon mode energy, 
the nearly free electron states start mixing with the ir 
bands. This leads to a gap of ±uj a around the Fermi level 
in the nearly free electron channel. Outside this gap, the 
nearly free electron DOS recovers the usual "V" -shape of 
graphene's DOS. (See Fig. 2, lower panel.) Indeed, the 
shape of the DOS in this channel is very similar to the 
gapped spectra found in STS on graphene. [2] 

Comparing the scales in the upper and lower part of 
Fig. 2 one sees, that the DOS in the nearly free elec- 
tron channel is for all coupling constants considered, here, 
smaller than the total DOS by a factor of less than 1/50. 
In STS, however, this factor of 1 /50 is by far overcompen- 
sated by the factor of order 10 4 coming from the different 
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Figure 2: (Color online) DOS for different coupling strengths 
A = 0, 0.3, 0.5, 0.7 eV. Upper panel: Total DOS. Lower 
panel: DOS of nearly free electron band. 
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Figure 3: (Color online) DOS of nearly free electron states 
for different chemical potentials fi — —0.4, —0.2, 0.0, 0.2 
and 0.4 eV and the coupling strength A = 0.5 eV fixed to the 
experimental value [2]. The curves are vertically displaced for 
clarity.. 



tunneling matrix elements in presence or absence of K/K' 
out of plane phonons. 

As the STM spectra are dominated by the nearly free 
electron channel, it is illustrative to study Nr as a func- 
tion of the chemical potential, as shown in Fig. 3. The 
gap of Lu a around the Fermi level is independent of dop- 
ing, while the minimum caused by nearly free electron 



states coupled with 7r-states at the Dirac point shifts with 
the chemical potential. Experimentally, the chemical po- 
tential is shifted with the gate voltage [2, 3]. So, the 
spectra calculated, here, and the experimental STS from 
Rcf. [2] exhibit the same characteristic behavior. 

The good agreement of the calculated tunneling DOS 
and the experimental STM spectra firstly proves the cor- 
rectness of the phonon mediated tunneling mechanism 
proposed in Ref. [2]. Moreover, it becomes clear that 
the phenomenon of strongly enhanced tunneling currents 
may occur in various contexts in graphene. No inelastic 
effects are required: The DFT calculations presented in 
the first part of this article have a static lattice and show 
the enhancement of the local DOS in vacuum in presence 
of frozen K/K' phonons. The only effect important for 
such an enhancement is the mixing of the Dirac electron 
states at the Brillouin zone corners K/K' with the nearly 
free electron states at T. Therefore, any out-of-plane cor- 
rugations on scale of the graphene lattice constant will 
increase the tunneling density of states significantly. This 
has important consequences for STM experiments: The 
"visibility" of impurity states at low energies [18] will be 
strongly enhanced as soon as the impurity causes out 
of plane distortions of the graphene lattice STM images 
below the phonon threshold[19] can be expected to look 
much more inhomogeneous than those taken at higher 
bias voltages. 
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